%--------------------------------------------------------------------------
% OPTION PRICING
%--------------------------------------------------------------------------

k_vec       = repmat(k',[nS,1]);
k_vec       = k_vec(:);

kb_vec      = repmat(k_bar,[nK,1]);
kb_vec      = kb_vec(:); % [nS*nK,1]

mu_vec      = repmat(mu_XQ,[nK,1]);
mu_vec      = mu_vec(:); % [nS*nK,1]
sig_vec     = repmat(sig_X,[nK,1]);
sig_vec     = sig_vec(:); % [nS*nK,1]
Q_mat       = repmat(Q,[nK,1]); % [nS*nK,nS]
Rf_mat      = repmat(Rf(:,1),[nK,1]);
Rf_mat      = Rf_mat(:);

default     = reshape( k' <= k_def ,[nK*nS,1]);
issuance    = reshape( k' >= k_iss ,[nK*nS,1]);

Put         = zeros(nS,nK);
PV_D        = zeros(nS,nK);
PV_SD       = zeros(nS,nK);

[Cheb_nodes,~] = qnwcheb(nE,-1,1);
a           = -5;
b           = 5;
Eshocks     = (b-a)*(Cheb_nodes'+1)/2 + a;
temp        = pi*(b-a)/(2*nE);
weights     = temp.*sqrt(1-Cheb_nodes'.^2);
probE       = normpdf(Eshocks).*weights;

% kappa transition
g               = exp(mu_vec + sig_vec.*Eshocks);
kp              = k_vec  .* g; %
kp_issuance     = kb_vec .* g; %
kp(issuance,:)  = kp_issuance(issuance,:);

for is=1:nS
    lamSx_fit   = griddedInterpolant(k,lamSx(is,:)','linear');
    lamSx_fine(is,:) = lamSx_fit(k)';
end

for is=1:nS % future state
    lamS_fit    = griddedInterpolant(k,lamS(is,:)','linear');
    lamSx_fit   = griddedInterpolant(k,lamSx(is,:)','linear');
    %Div_fit     = griddedInterpolant(k,Div(is,:)','linear');
    SDiv_fit    = griddedInterpolant(k,SDiv(is,:)','linear');
    
    
    lamSx_1     = lamSx_fit(kp);
    lamS_1      = lamS_fit(kp);
    %Div_1       = Div_fit(kp);
    SDiv_1       = SDiv_fit(kp);
    
    cash_flow   = max(0, lamSx_fine(:).*money(:) - g.*SDiv_1 - g.*lamSx_1 );
    
    Put         = Put + reshape(Q_mat(:,is).*sum(cash_flow.*probE,2)./Rf_mat,[nS,nK]); % [nS,nK] %%%
    
    Div_1       = lamS_1 - lamSx_1 - SDiv_1;
    
    PV_D        = PV_D  + reshape(Q_mat(:,is).*sum( g.*Div_1.*probE ,2)./Rf_mat ,[nS,nK]); % no dividend for issuance
end

Put(default) = NaN;


% invert option prices
opt     = optimset('TolX',1e-6,'Display','off'); % for solving BS for IV
BSIV    = NaN(nS,nK);

for ik=2:nK
    for is=1:nS
        if k(ik)>k_def(is)
            if Put(is,ik)~=0
                if ik>2
                    x0 = BSIV(is,ik-1);
                    if isnan(x0), x0 = 0.3; end
                else
                    x0 = 0.3;
                end
                BSIV(is,ik) = fzero(@(vol)BS_put(lamSx_fine(is,ik)-PV_D(is,ik),money(is,ik)*lamSx_fine(is,ik),log(Rf(is,1))*12,1/12,vol)-Put(is,ik),x0,opt);
            end
        end
    end
end
